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Abstract 

We  use  a  multiscale  natural  pixel  type  representation  of  an  object,  originally  developed  for 
incomplete  data  problems,  to  construct  nearly  orthonormal  basis  functions.  The  coefficients  of 
expansion  of  an  object  in  these  basis  functions  are  obtained  as  the  1-D  wavelet  transform  of  the 
(strip  integral)  projections  of  the  object.  This  enables  us  to  formulate  a  multiscale  tomographic 
reconstruction  technique  wherein  the  object  is  reconstructed  at  multiple  scales  or  resolutions. 

A  complete  reconstruction  is  obtained  by  combining  the  reconstructions  at  different  scales. 

The  nearly  orthonormal  behavior  of  the  basis  functions  results  in  a  system  matrix,  relating 
the  input  (the  object  coefficients)  and  the  output  (the  projection  data),  which  is  extremely 
sparse.  The  system  matrix,  in  addition  to  being  sparse,  is  well-conditioned  and  has  a  symmetric 
block- Toeplitz  structure  if  the  angular  projections  are  uniformly  spaced  between  0°  and  180°. 

Fast  inversion  algorithms  exist  for  these  matrices.  The  multiscale  reconstruction  technique  can 
find  applications  in  object  feature  recognition  directly  from  projection  data,  tackling  ill-posed 
imaging  problems  where  the  projection  data  are  incomplete  and/or  noisy,  and  construction  of 
multiscale  stochastic  models  for  which  fast  estimation  algorithms  exist.  In  this  paper,  we  include 
examples  illustrating  the  above  applications  of  our  multiscale  reconstruction  technique. 

’This  work  was  supported  by  the  Office  of  Naval  Research  under  grant  N00014-91-J-1004,  by  the  US  Army 
Research  Office  under  grant  DAAL03-92-G-0115,  and  by  the  Air  Force  Office  of  Scientific  Research  under  grant 
F49620-92-J-0002. 
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1  Introduction 


The  conventioneil,  and  most  commonly  used,  method  for  image  reconstruction  from  tomographic 
projections  is  a  transform  method  called  Convolution  Back-Projection  (CBP)  [1].  The  CBP  recon¬ 
struction,  though  fast,  is  not  suitable  for  imaging  problems  where  the  projection  data  are  incomplete 
(limited  angle  and/or  truncated  projections)  [2,  7]  or  noisy  and  where  the  fimdamental  interest  is 
not  in  the  actual  pixel  values  themselves,  but  rather  in  something  that  is  derived  from  these,  such 
as  averages,  boimdaries  [4]  etc.  These  problems  are  encotmtered  in  many  applications  in  medicine, 
non- destructive  testing,  oceanography  and  surveillance.  Here  we  present  a  transform  method  for 
image  reconstruction  where  we  work  in  a  multiscale  transform  space.  This  multiscale  reconstruc¬ 
tion  method  provides  a  framework  that  has  the  potential  of  overcoming  the  above  limitations  of 
the  CBP  reconstruction.  Specifically,  we  use  a  natural  pixel  type  representation  of  an  object  [5,  6] 
to  construct  nearly  orthonormal  basis  functions.  The  coefficients  of  expansion  of  the  object  in 
these  basis  functions  can  be  computed  from  the  projection  (i.e.  the  strip  integral)  data  by  using 
the  wavelet  transform. 

The  natural  pixel  method  [5,  6]  expands  the  object  in  the  same  basis  functions  along  which  the 
projection  data  are  collected,  thereby  using  a  basis  representation  for  the  object  that  is  closer  to  the 
measurement  domain  than  the  standard  rectangular  pixel  basis.  The  natural  pixel  representation 
results  in  a  matrix  based  reconstruction  method.  Since  matrix  based  reconstruction  methods  do 
not  utilize  the  space  invariance  assumption  of  the  CBP,  the  natural  pixel  reconstruction  is  devoid  of 
the  many  limited  data  artifacts  present  in  the  CBP  reconstruction.  The  disadvantage  of  the  natural 
pixel  reconstruction,  or  the  matrix  based  reconstruction  methods  in  general,  is  that  inversion  of 
very  large,  ill-conditioned  matrices  is  required. 

In  this  paper,  we  bmld  on  the  natural  pixel  approach  by  using  wavelet  bases  to  transform  the  natmal 
pixel  basis  functions.  The  use  of  the  wavelet  bases  enables  us  to  formulate  a  multiscale  tomographic 
reconstruction  technique  wherein  the  object  is  reconstructed  at  multiple  scales  or  resolutions.  The 
coarser  scales  contain  the  low  frequency  (i.e.  the  low  resolution)  information  about  the  reconstructed 
object  and  the  finer  scales  contain  the  high  frequency  (i.e.  the  high  resolution)  information.  The 
standard  reconstruction  is  obtained  by  combining  the  reconstructions  at  different  scales.  The 
natural  pixel  transformation  matrix,  relating  the  input  (the  object  coefficients)  and  the  output 
(the  projection  data),  is  full.  The  use  of  wavelet  bases,  in  addition  to  providing  a  multiscale 
framework,  restdts  in  this  transformation  matrix  being  sparse.  Also,  the  mxiltiscale  framework 
allows  us  to  use  simple  geometric  argmnents  to  partition  the  multiscale  transformation  matrix  such 
that  the  reconstruction  method  requires  only  the  inversion  of  a  sparse  and  well- conditioned  matrix 
which  is  symmetric  block-Toeplitz  if  the  projections  are  uniformly  spaced  between  0°  and  180°. 
Note  that  fast  inversion  algorithms  exist  for  these  matrices. 

Many  imaging  problems  are  iU-posed  in  the  sense  that  we  wish  to  reconstruct  more  degrees  of 
freedom  than  exist  in  the  data.  Noisy  and/or  incomplete  projections  (as  in  low  dose  medical  imag¬ 
ing,  oceanography,  and  in  several  applications  of  nondestructive  testing  of  materials)  make  the 
reconstruction  problem  ill-posed.  In  these  problems  the  lower  resolution  (i.e.  the  coarser  scale) 
reconstructions  are  more  reliable  than  their  higher  resolution  counterparts.  The  multiscale  re¬ 
construction  approach  provides  estimates  of  the  field  (or  the  object)  at  a  variety  of  resolutions 
thus  providing  a  natural  framework  for  explicitly  assessing  the  resolution-accuracy  tradeoff.  We 
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show,  through  examples,  that  for  an  ill-posed  reconstruction  problem,  some  regularization  can 
be  achieved  by  only  combining  coarser  scale  reconstructions  instead  of  reconstructions  at  all  the 
scales.  For  noisy  data  problems,  we  specialize  our  multiscale  reconstruction  method  to  yield  fast 
MAP  multiscale  reconstruction  estimates  corresponding  to  a  chosen  prior  on  the  multiscale  object 
coefficients. 

Apart  from  this,  the  multiscale  framework  is  useful  in  many  other  situations.  It  lends  itself  naturally 
to  situations  where  the  projection  data  are  collected  at  multiple  resolutions.  It  also  provides  the 
ability  to  reconstruct  different  parts  of  the  object  at  different  resolutions.  This  could  be  particularly 
useful,  for  example,  if  the  projection  data  is  gathered  with  spatially  varying  density  or  if  the  desire 
is  to  focus  computational  resources  at  certain  points  of  interest. 

Finally,  the  multiscale  reconstruction  method  provides  a  means  of  object  feature  extraction  directly 
from  the  projection  data.  For  instance,  if  we  are  only  interested  in  imaging  high  frequency  details 
within  the  object  (for  example,  botmdaries^),  then  we  could  directly  obtain  these  features  by  using 
only  the  finer  scale  information  in  the  data.  Similarly,  if  averages^  of  the  pixel  values  in  a  region  are 
desired,  then  only  the  coarser  scale  information  is  needed.  Using  conventional,  ad  hoc  techniques, 
we  would  first  have  to  reconstruct  the  object  and  then  use  post-processing  to  extract  such  features. 
The  disadvantages  of  the  conventional  approach  are  that  (a)  a  good  quality  reconstruction  (and 
hence  a  complete  projection  data  set)  is  required  for  post-processing  [4],  and  (b)  even  though  the 
noise  in  projection  measurements  is  white,  the  noise  in  the  reconstructed  image  is  colored,  and 
hence  the  necessity  of  using  a  2-D  whitening  filter  on  the  reconstructed  image  [3] . 

The  paper  is  organized  as  follows.  In  Section  2  we  describe  the  tomographic  reconstruction  problem. 
In  Section  3  we  describe  the  natural  pixel  reconstruction  technique  and  explore  it’s  relationship 
with  conventional  reconstruction  methods.  In  Section  4,  starting  from  the  natural  pixel  object 
representation,  we  develop  the  theory  behind  our  multiscale  reconstruction  method.  We  present 
some  sample  reconstructions  in  Section  5,  and  present  conclusions  in  Section  6.  Appendices  1-5 
contain  mathematical  details  of  our  multiscale  reconstruction  method,  and  Appendix  6  stunmarizes 
the  mathematical  notations  used  throughout  this  paper. 


2  The  Tomographic  Reconstruction  Problem 


In  non- diffraction  tomography,  the  goal  is  to  reconstruct  an  object  or  a  field,  /,  from  the  projection 
data  [1].  For  a  parallel-beam  imaging  geometry,  the  projection  data  consists  of  parallel,  non¬ 
overlapping  strip  integrals  through  the  object  at  various  angles  (refer  to  Figure  1).  Each  angular 
position  corresponds  to  a  specific  source-detector  orientation.  Thus  the  projection  data  can  be 
mathematically  represented  as  i/fe„,  where  k  —  1,...,M  (M  is  the  number  of  uniformly  spaced 
angular  positions  between  0  and  180°),  and  n  =  1, . . . ,  N  {N  is  the  ntunber  of  parallel  strip  integrals 
in  each  angular  position).  Furthermore,  we  call  the  indicator  function  of  the  strip  labeled  by 

’Boundary  detection  has  applications  in  medical  imaging  (detection  and  outlining  of  boundaries  of  organs  and 
tumors),  nondestructive  testing,  oceanography  and  plant  physiological  studies  [4]. 

^Computation  of  average  values  in  a  region  of  the  field  has  applications  in  functional  medical  imaging. 
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Figure  1:  The  projection  measurements  for  an  object,  /  (shaded),  at  two  different  angular  positions 
(fe  =  1  and  k  =  2  respectively).  The  number  of  parallel  strip  integrals  in  each  angular  projection, 
N,  is  8  in  this  case.  Three  basis  functions,  <t>n,4>is,<j>2S,  which  are  the  indicator  functions  of  the 
corresponding  strips,  are  also  shown. 


(fc,n)  so  that  has  value  one  within  that  strip  and  zero  otherwise.  Given  this  notation, 

ykn=jj  f{u,v)^kn{u,v)dudv  k  =  n  =  l,...,iV  (1) 

n 

where  (u,  u)  are  the  rectangular  coordinates  and  the  integration  is  carried  over  a  region  of  interest 

n. 


Due  to  practical  considerations,  we  have  to  work  with  a  discretized  version  of  Equation  1.  This  is 
y=[y^{l)  y^{2)  ...  y^{k)  ...  £^(M)]^  =  ^/  (2) 

where  y(fc)  is  the  projection  data  set  collected  at  angle  k,  and  is  composed  of  the  N  strip  integrals 
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as  shown  below  ^ 

y{k)  -  [  vki  yk2  •  •  •  ykN  ]  •  (3) 

The  vector  /  is  formed  by  first  discretizing  f{u,  v)  on  a  rectangular  lattice  (we  will  assume  apxp 
square  pixel  lattice)  and  then  stacking  the  discretized  values  in  a  lexicographic  order  into  a  vector 

/=  [  A  A  /p2  (4) 

and  the  MN  X  p^  matrix  <!>  is  formed  by  first  stacking  aU  {<f>kn{u,  v),  ^  =  1, . . . ,  M;  n  =  1, . . . ,  N} 
into  a  vector  followed  by  discretization  in  u  and  v  (on  the  same  square  pixel  lattice  used  for  /)  and 
lexicographic  ordering  (to  match  that  in  /) 


(5) 


■  &  ' 
& 

<f>ll,2 

■  •  • 

(j>  = 

— 

4>12,1 

<l>\2,2 

-&N. 

<I>MN,2 

•  •  •  4>MN,p^ 

The  tomographic  reconstruction  problem  then  reduces  to  finding  an  estimate,  /,  of  the  discretized 
object,  /,  given  the  projection  data  y. 


3  The  Conventional  Reconstruction  Techniques 


In  this  section  we  discuss  two  conventional  reconstruction  techniques,  the  widely  used  convolu¬ 
tion  back-projection  (CBP)  reconstruction  technique,  and  the  natirral  pixel  (NP)  reconstruction 
technique  used  by  us  as  a  starting  point  for  our  multiscale  reconstruction. 


In  both  the  CBP  and  the  NP  reconstructions,  the  object  is  represented  as  a  linear  combination  of 
(non-orthogonal)  basis  functions  <j}  along  which  the  projection  data  are  collected.  Mathematically, 

M  N 

f{u,v)  ^kn<f>kn{u,v).  (6) 

k=l  n=l 


The  discretized  version  of  the  above  equation  is 

f  =  <Fx 


(7) 


where 

£=  [ 

In  the  above  equation,  x(k)  is 
coefficients  as  shown  below 


5.^(1)  ^^(2)  •  ■  •  ^(k)  •  •  ■  x^{M)  j  .  (8) 

the  object  coefficient  set  at  angle  k,  and  is  composed  of  the  N  object 


x{k)  = 


Xkl  Xk2 


(9) 


In  the  NP  reconstruction,  the  coefficients  of  expansion,  are  called  the  natural  pixels.  Equa¬ 
tion  7  can  also  be  interpreted  as  the  back-projection  operation  where  the  coefficients  Xkn  are 
back-projected  along  the  basis  functions  Thus  if  the  natural  pixels,  Xkn^  are  the  same  as  the 
ramp-filtered  projection  coefficients  of  the  CBP,  the  natural  pixel  and  the  CBP  reconstructions  will 
give  identical  results.  We  wiU  come  back  to  this  later. 
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3.1  The  convolution  back-projection  (CBP)  reconstruction 

The  CBP  reconstruction  assumes  a  space-invariant  point  spread  function  [7,  8,  9]  and  hence  is 
only  valid  for  complete  data  problems.  We  borrow  the  definition  of  a  complete  data  set  from 
Llacer.  According  to  Llacer  [7],  “a  complete  data  set  could  be  described  as  sufficient  number 
of  line  projections  at  a  sufficient  number  of  angtdar  increments  such  that  enough  independent 
measurements  are  made  to  allow  the  image  reconstruction  of  a  complete  bound  region.” 

In  the  CBP  reconstruction,  the  space  invariance  property  is  utilized  to  calculate  x{k)  at  each 
angular  position,  k,  by  convolving  the  projection  data  at  that  particul^lr  angular  position  with  the 
inverse  Fourier  treuisform  of  a  ramp  filter  [1].  Thus,  for  a  fixed  angular  position  k, 

x{k)  =  F-'^{p)*y{k) 

=  TaV^k)  (10) 

where  p  is  a  vector  containing  the  1-D  ramp  filter  coefficients,  F~^{p)  represents  the  1-D  inverse 
fourier  transform  of  the  1-D  ramp  filter,  *  refers  to  1-D  convolution  in  the  spatial  domain,  and  the 
matrix  performs  the  operation  F~^{p)*.  Let 

T  =  lMxM®Ta  (11) 

where  ImxM  is  an  M  X  M  identity  matrix  and  (gi  refers  to  the  Kronecker  product.  The  above 
equation  implies  that  T  is  a  block-diagonal  matrix  with  M  blocks  along  the  diagonal,  all  equal  to 
Tq.  Thus,  for  the  CBP  reconstruction,  the  necessary  steps  are  as  follows^  (a)  computation  of  the 
basis  coefficients  x  from  projection  data  y  by  the  equation 

x-Ty_  (12) 

where  T  is  as  defined  in  Equation  11,  and  (b)  back-projection  according  to 

f_-(l)'^x  =  (j)^Ty.  (13) 

Note  that  if  all  elements  of  p  are  replaced  by  unity,  then  Equation  13  reduces  to  simple  back- 
projection  reconstruction  (without  any  filtering).  This  reconstruction  is  approximately  equal  to  the 
true  object,  f{r,d)  in  radial  coordinates,  blurred  with  a  (circularly  symmetric)  l/r  point  spread 
function. 


3.2  The  natural  pixel  (NP)  reconstruction 

In  the  NP  reconstruction  [5,  6],  the  space- invariance  assumption  of  the  point  spread  function  is 
not  used  to  find  the  natural  pixels,  from  the  projection  data,  y,  and  hence  it’s  applicability 
to  incomplete  data  problems.  Equations  2  and  7  are  used  instead.  Thus  the  CBP  reconstruction 
can  be  viewed  as  a  special  case  of  the  NP  reconstruction  under  space-invariance  or  complete  data 
assumptions.  By  substituting  f  from  Equation  7  for  /  in  Equation  2,  we  get 

y  =  =  Gx  where  G  =  (jxjP^ .  (14) 

®In  practice,  for  regularization  purposes,  a  roll-off  is  applied  to  the  ramp  filter.  This  operation  would  include  an 
additional  matrix  Q  and  so  F~^{p)  will  be  replaced  by  F~^(Qp)  [1]. 
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Figure  2:  The  elements  of  the  matrix  G  are  the  areas  of  intersection  of  various  strips.  One  such 
area  of  intersection,  corresponding  to  two  strips  delineated  by  bold  lines,  is  shown  shaded.  The 
matrix  G  is  full  as  most  of  these  areas  are  non-zero. 


The  elements  of  the  matrix  G  are  the  areas  of  intersection  of  the  strips  defined  by  the  basis  functions 
4>.  Most  of  these  areas  are  not  zero  and  hence  the  matrix  G  is  full  (refer  to  Figmre  2).  It  can  be 
easily  shown  that  G  is  symmetric  block- Toeplitz  if  the  M  angular  projections  are  uniformly  spaced 
between  0°  and  180°,  even  though  G  may  be  quite  ill-conditioned  in  general,  as  we  will  see  later. 
Figure  7,  top,  shows  G  for  an  imaging  geometry  with  32  angular  projections  (i.e.  M  =  32)  and 
with  32  strip  integrals  (i.e.  N  =  32)  in  each  angular  projection. 

The  advantage  of  the  NP  reconstruction  over  the  CBP  reconstruction  is  that  (a)  since  the  matrix 
G  can  be  calculated  or  measured  for  each  specific  geometry,  the  reconstruction  can  be  customized 
for  any  imaging  system,  and  (b)  the  point  spread  fimction  does  not  have  to  be  space  invariant 
and  hence  a  complete  set  of  angular  projection  data  is  not  reqxured  as  in  CBP.  However,  the  NP 
reconstruction  suffers  from  some  disadvantages  as  well,  namely  that  (a)  the  large  size  of  the  matrix 
G  {MN  X  MN)  imposes  a  limitation  on  storage  and  speed,  and  (b)  the  matrix  G  can  be  quite 
ill-conditioned  for  some  imaging  geometries  (see  the  following  discussion  and  [7,  8,  9]). 

In  the  next  section,  we  use  wavelet  bases  to  transform  the  NP  basis  functions,  (f>,  into  the  multiscale 
framework.  As  described  in  the  Introduction,  we  gain  many  important  featnres  in  going  to  the 
multiscale  framework.  In  addition,  we  overcome  the  above  limitations  of  the  NP  reconstruction. 
The  multiscale  system  matrix,  Pg,  is  related  to  (?  by  a  similarity  transformation  and  is  sparse.  This 
and  the  fact  that  Pg  is  block- Toeplitz  for  the  case  of  uniformly  spaced  angvdar  projections,  offsets  the 
first  disadvantage  of  the  NP  reconstruction.  Also,  the  geometry  of  the  imaging  system  is  naturally 
captured  in  P^.  We  exploit  this  featme  to  partition  Pg  into  weU  and  ill-conditioned  matrices  such 
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that  our  reconstruction  procedure  requires  inversion  of  only  well- conditioned  matrices.  This  offsets 
the  second  disadvantage  of  the  NP  reconstruction. 

Connection  with  the  least-squares  reconstruction 


If  the  matrix  has  full  row  rank,  then  the  NP  matrix  G  is  invertible  and  in  this  case  the  NP 
solution  is  given  by  /  =  This  is  just  a  special  case  of  the  minimum  2-norm  least  squares 

solution  to  y  =  <f>f  (Equation  2)  and  given  by  [12] 

Ils  =  (15) 

where  refers  to  the  pseudo-inverse  of  (j).  However,  as  we  discuss  next,  it  is  not  possible  in  general 
to  invert  G  even  when  <j)  has  full  row  rank.  This  is  because  the  matrix  G  is  ill-conditioned  due  to 
the  inherent  non-uniqueness  in  the  NP  representation. 

Non-uniqueness  in  the  NP  representation 


Because  the  NP  object  representation  is  tied  to  the  collection  of  projection  data,  there  is  an  in¬ 
herent  non-uniqueness  in  the  representation  which  results  in  G  not  having  full  rank  or  being  badly 
conditioned.  Specifically,  G  =  has  full  rank  (i.e.  is  invertible)  if  and  only  if  <f)  has  full  row 
rank.  Recall  from  Equation  5  that  the  rows  of  4>  aje  composed  of  discretized  basis  functions, 
•  •  •  ’  ^n’  ■  ■  ■  ’  &IN'  ^  rank  if  none  of  these  basis  functions  can  be  written  as  a 

linear  combination  of  the  others.  This  is  not  true  if  two  disjoint  subsets  of  span  the  same 

object  subspace.  As  an  example,  assume  that  all  angidar  projection  sets  exactly  span  the  entire 
object  space.  Then  obviously 

N  N 

-  Y,  ik^n  =  ^  V  fci  7^  and  *1,^2  =  (16) 

n=l  n=l 

because  the  DC  term  in  the  object  appears  in  all  M  angular  projections.  Thus,  in  this  example,  <f> 
does  not  have  full  row  rank,  resulting  in  a  singular  G.  Throughout  this  paper,  we  assume  that  p  =  N 
(i.e.  the  object  is  discretized  mapxp  =  NxN  rectangular  pixel  lattice,  where  N  =  number  of 
strip  integrals  in  each  angular  projection)  and  that  the  M  angular  projections  are  uniformly  spaced 
between  0  and  180°.  Thus  <f>  and  G  are  MN  x  MN  square  matrices.  In  this  geometry.  Equation  16 
is  true  only  if  ki  and  k2  correspond  to  projections  at  0°  and  90°,  since  only  projections  at  these 
angles  exactly  span  the  entire  object  space.  Now,  angtdar  projections  at  0°  and  90°  exist  if  and 
only  if  M  is  even  because  of  our  assumption  of  uniformly  spaced  angular  projections.  Thus  if  M 
is  even,  (f)  has  rank  MN  —  1,  otherwise  it  has  a  full  rank  of  MN.  Similarly,  it  can  be  shown  that 
G  drops  a  rank  for  each  pair  of  projections  90°  apart.  There  are  projection  pairs  90°  apart  if  and 
only  if  M  is  even,  in  which  case  there  are  exactly  M/2  of  them.  Thus  if  M  is  even,  G  has  a  rank 
of  MN  —  M/2,  otherwise  it  has  a  full  rank  of  MN. 

Thus,  if  M  is  even,  it  is  not  possible  to  uniquely  solve  for  ®  in  y  =  Gx.  Because  of  numerical 
issues  (we  use  7  digit  precision  for  calculation  of  elements  of  G),  the  matrix  G  has  full  rank  for  even 
as  well  as  odd  M.  The  condition  number  is  a  good  measure  to  use  in  this  case.  In  Figure  3,  we 
show  the  condition  numbers  of  G  for  imaging  geometry  with  M  =  N  and  for  values  of  M  ranging 
from  8  to  18.  Notice  the  steep  fluctuation  in  condition  numbers  between  odd  and  even  values  of 
M,  and  the  steady  increase  of  the  condition  number  for  increasingly  odd  values  of  M.  This  is 
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because  projection  pairs  approach  being  90°  apart  as  M  is  increased  in  this  case.  Thus,  due  to  this 
ill- conditioning,  it  is  not  possible  to  uniquely  solve  for  xmy_  =  Gx  even  for  odd  values  of  M.  We 
show,  in  Section  4,  that  this  non-uniqueness  problem  of  the  NP  representation  is  easily  dealt  with 
in  the  multiscale  framework.  As  a  restdt,  the  system  matrix  in  the  multiscale  framework  can  be 
partitioned  into  well  and  ill-conditioned  matrices  such  that  the  reconstruction  requires  inversion  of 
only  well- conditioned  matrices. 

The  dropping  of  the  rank  of  G  for  each  projection  angle  pair  90°  apart  is  illustrated  in  Figure  4. 
For  simpHcity,  consider  a  2  x  2  object 

l=[fi  fz  fz  Uf  (17) 

and  projection  measurements 

y  =  [s/ll  S/12  S/21  S/22  ]  (18) 

[  /l  +  /s  fz  +  fi  fl  +  fz  /s  +  /4  j  •  (19) 
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yii 

yi2 

f2 

^21 

f3 

f4 

^22 

Figure  4:  Imaging  geometry  where  the  object,  /,  is  discretized  on  a  =  4  pixel  lattice.  Here 
M  —  N  —  2  with  the  projections  being  90°  apart. 


Here 


has  a  rank  of  3,  and 


also  has  a  rank  of  3. 


10  10 
0  10  1 
110  0 
0  0  11 


G  =  (fxjp- 


2  0  11 
0  2  11 
112  0 
110  2 


(20) 


(21) 


4  The  Multiscale  Reconstruction 


The  multiscale  reconstruction  is  motivated  by  the  following  observations.  In  the  previous  section 
we  saw  that  the  NP  reconstruction  involves  solving  for  the  natural  pixels,  z,  from  the  projection 
data,  y,  which  are  related  by  y  =  Gx  and  then  back-projecting  x,  £=  4)^x.  Also,  the  elements 
of  matrix  G  are  the  areas  of  intersection  of  the  various  strips  along  which  the  projection  data  are 
collected.  Suppose  now  that  the  strips  are  of  the  form  shown  in  Figure  5.  Each  strip  is  a  linear 
combination  of  two  NP  strips,  one  given  a  positive  weight  and  the  other  negative.  The  new  matrix 
relating  x  and  y,  according  to  the  above  choice  of  strips,  will  have  as  its  elements  the  areas  of 
intersections  of  the  newly  defined  strips.  It  is  clear  from  Figure  5  that  most  of  these  elements  will 
be  zero  due  to  the  cancellation  of  the  positive  and  the  negative  terms.  Only  those  elements  that 
correspond  to  strip  intersections  near  the  edge  of  the  field-of-view  wiU  be  non-zero.  Thus  one  can 
expect  this  new  matrix  to  be  sparse  with  the  degree  of  sparsity  increasing  with  the  size  of  the 
field-of-view. 
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+ 


>(• 


Figure  5:  Each  of  the  two  delineated  strips  from  figure  2  are  broken  down  into  two  substrips,  having 
a  positive  and  a  negative  weight  respectively.  The  area  of  intersection  of  the  two  strips  is  zero  in 
this  case  due  to  the  cancellation  of  the  positive  (lightly  shaded)  and  the  negative  (darkly  shaded) 
terms. 


The  above  redefinition  of  the  strips  with  positive  and  negative  weights  is  reminiscent  of  the  finest 
scale  contribution  from  the  Haar  transform.  However,  cin  important  point  to  note  here  is  that  the 
Haar  transform  is  taken  only  in  one  direction,  i.e.  the  direction  perpendicular  to  the  long  axis  of 
the  strip.  As  an  example,  for  a  projection  at  a  fixed  angle,  A:  =  1,  and  consisting  of  eight  strips  (i.e. 
N  =  8),  the  full  Haar  transform  wiU  look  as  shown  in  Figure  6.  A  notion  of  scale  emerges  from  the 
use  of  the  Haar  transform.  The  original  strips  have  been  broken  down  into  strips  at  multiple  scales 
having  positive  and  negative  weights.  The  finest  scale  involves  strips  that  have  twice  the  width 
of  the  original  strips  and  the  coarsest  scale  involves  strips  that  have  eight  times  the  width  of  the 
original  strips.  We  wiU  call  the  above  transformation  of  the  strips  the  natural  wavelet  transform 
because  of  the  adaptation  to  the  natural  pixel  representation. 

Let  the  new  strips  be  defined  by  the  basis  functions  (recall  that  the  original  strips  were 

defined  by  basis  functions  {^^}).  We  use  lowercase  for  the  original  NP  basis  functions  and  up¬ 
percase  for  the  multiscale  basis  functions.  Note  that  each  multiscale  basis  function  is  labeled 
by  three  subscripts,  ksn,  where  k  indicates  the  angle  of  the  projection  {k  =  1, . . . ,  M),  s  indicates 
the  scale  (s  =  1,2, ...,1  -f  hiAr/ln2;  N  —  2*,  i  £  Z"^),  and  n  indicates  the  shift  within  the  scale 
(n  =  1, 2, . . . ,  n,;  n,  =  iV/2*  for  s  =  1, 2, . . .,  In  AT/ In  2  and  n,  =  1  for  s  =  1  -|-  In  A^/ln2).  In  the 
future,  we  will  refer  to  the  finest  scale  (s  =  1)  as  h  (for  high),  the  coarsest  scale  (s  =  In  A^/ln2)  as 
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The  natural  pixel  basis  functions  corresponding  to  projection  k=l 


Scale  1  decomposition 


Scale  3  decomposition 


<I> 


131 


Scale  2  decomposition 
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141 


Figure  6:  The  basis  functions  (or  equivalently  the  strips  defined  by  them)  are  decomposed 

into  functions  at  multiple  scales  given  by  {§i5n}- 


I  (for  low),  and  the  DC  component  (s  =  1  +  luiV/  ln2)  as  d  (for  dc).  Thus  s  =  h,h+  I, . .  .,l,d. 

Let  Wa  be  a  matrix  representation  of  the  linear  operator  that  performs  a  1-D  orthogonal  wavelet 
decomposition  on  a  discrete  sequence^.  Then,  at  a  fixed  angular  position  k,  our  new  multiscale 

^We  use  the  Discrete  Periodic  Wavelet  Transform  (DPWT)  [11]  with  the  Haar  and  the  various  Daubechies  com¬ 
pactly  supported  wavelets  [10]  for  our  reconstructions.  A  Daubechies  compactly  supported  wavelet  Dm  is  character¬ 
ized  by  m  vanishing  moments. 
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basis  functions  and  the  previous  NP  basis  functions  are  related  by 


■  ILl  ■ 

—khni. 

II 

; 

Tiir 

: 

.  ^idi  . 

(22) 


with  W~^  =  NW^  because  of  our  choice  of  orthogonal  wavelet  basis  normalized  such  that  the 
absolute  values  of  elements  in  the  multiscale  counterpart  of  matrix  G  lie  between  0  and  1  irrespective 
of  the  size  of  the  imaging  system.  As  an  example,  for  AT  =  8  and  using  the  Haar  basis  (Figure  6), 
the  bases  at  a  single  angle  are  related  by 


where 


Wa 


L  :^41  J 


% 
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(23) 


(24) 


Let 

W  =  ImxM  <8»  Wa  (25) 

i.e.  IF  is  a  block- diagonal  matrix  with  M  blocks  along  the  diagonal,  all  equal  to  Wa-  Then 

$  =  W(f>  (26) 

with  W~^  =  NW^,  where 

$  =  [  ^i;ii  •  •  •  .  .  ^Mhi  •  •  •  l-Mdi  ]  (27) 
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and  (j)  is  as  defined  in  Equation  5. 

Recall  the  natural  pixel  equations  (2),  (7)  and  (14).  The  multiscale  decomposition  of  vectors  x  and 
y,  given  by  ^  and  ^  respectively,  is  defined  as 


=  NWx 


(28) 


and 


»  (29) 

where  we  use  the  subscript  6  to  indicate  ordering  according  to  projection  angles  corresponding  to 
our  choice  of  W  in  Equation  25.  Thus,  to  illustrate  9  ordering. 


xb 


[&)  iSm  ...  i^(k)  ...  £(M)f 


(30) 


where  is  the  1-D  wavelet  transform  of  the  projection  data  set  at  angle  k,  yi^k),  and  is  composed 

of  the  elements  at  different  scales  as  shown  below 

t^{k)  ^  Way{k)  -  jJ^{k,h)  ^{k,h+l)  ...  i^{k,j)  ...  J^{k,d)  .  (31) 

In  the  above,  ^(fc,  j)  is  the  vector  containing  the  jth.  scale  components  of  the  1-D  wavelet  transform 
of  y{k) 

7 

3^(A:,j)  =  [  tpkji  •••  ] 


with 


^  for  j^d 
1  for  j  =  d. 


(32) 

(33) 

(34) 


Similarly 

£,  =  [«J(i)  £j(2)  SW  IJW 

where  ^{k)  is  the  1-D  wavelet  transform  of  the  object  coefficient  set  at  angle  k,  x{k),  and  is 
composed  of  the  elements  at  different  scales  as  shown  below 

iT 


Uk)  =  NW,x{k)^  ^ik,h)  ii{k,h+l) 


(35) 


In  the  above,  ^{k,  j)  is  the  vector  containing  the  jth  scale  components  of  the  1-D  wavelet  transform 
of  x{k) 

-  ,  T* 

(36) 


^{k^j)  —  ^kjl  ^kj2  •••  ^kjnj  j 


with 


N 

.  =  i  V 


n,-  = 


for  j  ^  d 
1  for  j  =  d. 


(37) 


The  three  subscripts  fcsn,  in  'tpksn  ^.nd  ^ksn  above,  convey  the  same  meaning  and  take  on  the  same 
values  as  those  for 


In  the  multiscale  coordinates,  the  input-output  equation  (14)  changes  to 

y  =  Gx 

Wy  =  WGW-'^Wx 

^  =  {WGW^)i^ 

=  (38) 
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where  the  symmetric  block- Toeplitz  matrix  Tg  is  given  by 

Tg  =  WGW^  =  (39) 

The  above  equation  implies  that  the  elements  of  the  matrix  are  the  areas  of  intersection  of  the 
various  multiscale  basis  functions  We  know  from  our  earlier  discussion  that  F#  is  sparse.  The 
object  representation  equation  (7)  changes  to 

/  =  (f>Tx  =  4>'^W-'^Wx  =  =  iW(l>f^ 

=  (40) 


In  some  cases  it  is  more  intuitive  to  order  according  to  scales  rather  than  projection  angles.  We 
will  denote  the  scale  ordering  by  subscript  s.  Thus 


i  - 

(41) 

(42) 

(43) 

and 

F,  =  PVgP- 

=  PFgP^ 

(44) 

where  P  is  the  (orthogon^d)  permutation  matrix  which  changes  projection  angle  ordering  into  scale 
ordering.  Thus,  to  illustrate  s  ordering. 

3^  =  [ 

+  ••• 

(45) 

where  V’,(j)  contains  projection  data  at  scede  j  from  midtiscale  projection  data  sets  at  aU  angles 
and  is  given  by 

±,{3)  = 

1^(1, i)  ^{2,j)  .. 

•  ...  5^(M,i)]^. 

(46) 

Similarly 

'  i^{h)  i^{h  +  i)  ... 

if(i)  ...  Cil)  Cid)Y 

(47) 

where  contains  object  coefficients  at  sccde  j  from  object  coefficient  sets  at  aU  angles  and  is 

given  by 


1,(7)  =  [  •••  S{k,j) 


(48) 


As  explained  in  the  beginning  of  this  section,  we  expect  the  multiscale  matrices  to  be  sparse.  The 
number  of  non-zero  elements  in  F,  and  Fg  are  the  same  and  hence  the  degree  of  sparsity  achieved 
in  the  multiscale  framework  can  be  studied  using  either  of  these.  Figure  7,  bottom,  shows  F, 
for  an  imaging  geometry  with  32  angtdar  projections  (i.e.  M  =  32)  and  32  strip  integrals  in  each 
angular  projection  (i.e.  N  =  32).  We  have  used  the  Haar  wavelet  for  this  miiltiscale  representation. 
Comparing  with  Figure  7,  top,  which  shows  the  corresponding  NP  system  matrix,  G,  we  see  that  F, 
is  sparser  than  G.  From  the  figure,  most  of  the  non-zero  terms  in  F,  correspond  to  the  coarser  scale 
terms  where  the  edge  effects  are  more  pronounced.  It  was  claimed  in  the  previous  section  that  F, 
becomes  more  sparse  as  the  size  of  the  field-of-view  increases.  This  claim  is  validated  by  Figure  8 
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Figure  7:  Top:  Natural  pixel  system  matrix,  G\  Bottom:  Multiscale  natural  pixel  system  matrix, 
r,,  for  an  imaging  geometry  with  32  angular  projections  (i.e.  M  =  32)  and  with  32  strip  integrals 
(i.e.  N  —  32)  in  each  angular  projection.  The  Haar  wavelet  is  used  for  multiscaie  decomposition. 


where  we  plot  the  degree  of  sparsity  of  F,  as  a  function  of  the  parameter  M  {—  N).  We  measure  the 
degree  of  sparsity  by  the  percentage  of  elements  in  matrix  F,  which  are  below  a  certain  threshold 
(the  maximum  value  in  F,  is  normalized  to  1).  Figure  8  reports  the  sparsity  calculations  for  three 
different  values  of  threshold,  namely  0,  0.005  and  0.02.  It  is  empirically  observed  that  setting  all 
values  in  F,  below  a  threshold  of  0.02  to  zero,  makes  no  visible  difference  to  the  reconstructions. 
From  Figure  8,  we  see  that  for  the  case  of  M  =  N  —  128  and  a  threshold  of  0.02,  F,  is  98.75% 
sparse  (or,  equivalently,  1.25%  full).  In  Figure  9,  we  show  the  degree  of  sparsity  of  F^  or  F,,  for 
M  —  N  —  32,  achieved  by  the  Haar  wavelet,  and  the  Daubechies  wavelets  D3  and  Dg.  From  the 
figure  we  see  that  the  number  of  elements  that  are  exactly  zero  decrease  as  wavelets  with  larger 
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Figure  8:  The  degree  of  sparsity  (percentage  of  elements  below  a  threshold)  of  multiscale  system 
matrix,  T,,  as  a  function  of  M  {=N).  The  Haar  wavelet  is  used  here  for  multiscale  decomposition. 
The  majdmum  value  in  T,  is  normalized  to  1.  Setting  aU  elements  in  T,  below  0.02  to  zero  makes 
no  visible  difference  to  the  reconstructions. 


support  are  used.  But,  since  we  can  always  threshold  aU  elements  of  the  matrices  below  0.005  to 
zero  without  affecting  the  reconstructions,  the  effective  sparsity  achieved  by  the  Haar,  and  Ds 
is  approximately  the  same. 


Next,  we  tackle  the  non-uniqueness  issue  of  the  NP  representation.  We  had  seen  earlier  that,  with 
our  assmnption  of  p  =  iV  (i.e.  the  object  is  discretized  ma,pxp  =  NxN  rectangular  pixel 
lattice,  where  N  —  number  of  strip  integrals  in  each  angulcir  projection)  and  M  angular  projections 
uniformly  spaced  between  0  and  180°,  G  has  a  rank  of  MN  —  M/2  for  M  even  and  a  (full)  rank 
oi  MN  for  M  odd.  Tg  and  F,  have  the  same  rank  as  G.  Thus,  if  M  is  even,  it  is  not  possible  to 
uniquely  solve  for  ^(|^)  in  =  re|^(^  =  r,|^).  We  again  consider  the  same  2x2  example  as 
before  (Figure  4).  For  this  example,  using  the  Haar  basis. 
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Figure  9:  The  degree  of  sparsity  (percentage  of  elements  below  a  threshold)  of  multiscale  system 
matrix,  F,,  as  a  function  of  threshold,  for  different  wavelets.  Here  M  =  N  =  32.  The  maximum 
value  in  T,  is  normalized  to  1.  Setting  all  elements  in  F,  below  0.005  to  zero  maikes  no  visible 
difference  to  the  reconstructions. 


has  a  rank  of  3  because  the  two  DC  rows  (and  columns)  corresponding  to  the  two  projections  90° 
apart  are  identical.  In  general  for  M  projections,  where  M  is  even,  M/2  of  these  rows  and  columns 
wiU  be  identical  resulting  in  a  MN  -  M/2  rank  matrix.  This  suggests  a  recipe  for  partitioning  F, 
into  ill  and  well-conditioned  matrices.  Taking  advantage  of  this,  we  devise  an  approximation  to  the 
exact  equation  ^  =  F,^^  which  requires  inversion  of  a  relatively  weU-conditioned,  full  rank  matrix. 
This  approximate  reconstruction  is  vsdid  for  all  M,  even  or  odd,  and  yields  reconstructions  which 
are  almost  identical  to  the  CBP  reconstruction  for  the  complete  data  problem,  and  much  better 
for  the  incomplete  data  case.  We  partition  ^  according  to 


t’{d) 


‘  F.i 

r,2 

r,<i 

L- 

Lid) 


(50) 


where  F,i  is  an  M{N  —  1)  X  M{N  —  1)  symmetric  matrix,  F,2  is  an  M  x  M(N  —  1)  matrix,  and  Fjd 
is  an  M  X  M  symmetric  circuleint  matrix.  The  vectors  ^  (d)  and  ^^(d)  have  length  M  and  contain 
the  DC  terms.  The  vectors  and  have  length  M{N  —  1)  and  are  given  by 


3^-^  •••  t^(i)  •••  :^(0  (51) 
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and 

£^('•  +  1)  •••  £f«)  (52) 

The  —  in  the  subscript  for  and  indicates  that  these  vectors  do  not  contain  the  DC  terms. 

As  a  consequence  of  using  the  DPWT,  the  DC  (i.e.  the  last)  row  of  matrices  Wa  and  the  DC  basis 
functions  identical  irrespective  of  the  wavelet  used.  Thus,  if  we  assume  that  the  object  is 

fully  covered  by  aU  M  projections,  all  M  elements  of  ^(d),  which  are  the  total  number  of  counts 
in  the  corresponding  angular  projections,  are  equal  irrespective  of  the  wavelet  used.  Let  the  total 
number  of  cotmts  be  /x.  Thus 

(53) 

where  refers  to  a  vector  of  length  M  with  all  elements  equal  to  unity. 


As  mentioned  before,  T,  is  singular  for  even  M.  For  odd  M,  F,  is  invertible  even  though  it  may  be 
quite  ill-conditioned  because  of  projection  padrs  nearly  90°  apart.  In  both  cases,  M  being  even  or 
odd,  Tj!  is  well- conditioned.  For  even  M,  to  overcome  the  problem  of  non-invertibility  of  F,,  we 
assume  all  M  elements  of  the  vector  to  be  equal.  Let  this  value  be  c.  Then 


Thus 


’  t.-  ' 

'  r,i 

i,- 

.  n'^  . 

. 

cIm  . 

'  F,i 

F,2 

u 

(54) 


(55) 

(56) 


where  u  is  a  vector  of  length  M{N  -  1)  the  elements  of  which  are  the  row  sums  of  F^25  u  is 
a  vector  of  length  M  containing  the  row  sums  of  F,d.  Now,  the  elements  of  matrix  T,d  are  the 
areas  of  intersection  of  basis  function  {Jljfedi)  k  =  1,...,M},  and  because  of  our  assumption  of 
uniformly  spaced  angles,  F^d  is  a  circulant  matrix  (refer  to  Appendix  1).  Thus  the  row  sums  of 
F,d,  and  hence  all  M  elements  of  vector  u,  will  be  equal.  Let  this  value  be  a.  Thus 


U  =  0.1]^. 


(57) 


Now,  by  manipulation  of  Equations  50  -  57  (refer  to  Appendix  2  for  details),  our  multiscale 
reconstruction  method  consists  of  the  following  steps. 


1.  Form  the  matrix 


r.= 


■  r,i 

F,2 

F,d 

=  PVeP^'  =  P{WGW^)P'^ 


which  is  sparse. 

2.  Treinsform  the  projection  data  in  the  multiscale  basis 


= 


tJ-d) 


=  PWy. 


(58) 


(59) 
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3.  Calculate  the  multiscale  natural  pixels 


L-  =  ^  (60) 

where  the  constants  A  and  b  depend  only  on  the  parameters  of  the  imaging  system  and  the 
total  counts  in  the  object,  and  are  given  by 


—  aM  ‘ 


(61) 


and 


‘  =  (^K 


»1 


N 


-i. 
■ii  ^*1 2 


aM 


)■ 


(62) 


Note  that  Equations  61  and  62  require  the  inverse  of  Eji  which  is  a  full  rank,  relatively  well- 
conditioned  matrix®.  The  condition  number  of  T,!,  using  the  Haar  wavelet,  for  M  =  iV  =  16 
is  209  and  for  M  =  iV  =  32,  1288.  The  condition  number  of  F,,  on  the  other  hand,  is  3  x  10^ 
and  1.4  X  10®  respectively,  for  these  two  configurations.  Also,  since  elements  of  v  are  mostly 
negligible  (refer  to  Appendix  3),  a  good  approximation  to  Equations  61  -  62  is 


and 

and  thus 

4.  Back-project 

where 


A  ~ 


6  w  0 


L-  ^ 


f  =  +  Ji-i 

L  ~  7V2-' 


N't 


LiV2 


=  P$. 


(63) 

(64) 

(65) 

(66) 

(67) 


To  justify  the  above  reconstruction  algorithm  for  odd  M,  we  have  to  simply  justify  the  assumption 
that  aJl  M  elements  of  i^{d)  are  equal  for  odd  M.  This  is  done  in  Appendix  4. 

To  recapitulate,  the  miiltiscale  reconstruction  procedure  as  outlined  above  is  based  on  the  following 
two  assumptions. 


1.  The  elements  of  the  matrix  r,2  are  mostly  negligible.  This  is  empirically  seen  to  be  true  for 
the  Haar  and  the  compactly  supported  Daubechies  wavelets.  In  Appendix  3,  we  demonstrate 
the  validity  of  this  assumption  in  the  Haar  case  by  calculating  numerical  boimds  on  the 
absolute  values  of  the  elements  of  rj2. 

2.  AU  M  (DC)  elements  of  are  equal  if  the  object  is  entirely  covered  by  all  M  projections. 
We  show,  in  Appendix  4,  that  this  assumption  is  exactly  true  when  F,  has  full  rank,  i.e.  when 
M  is  odd,  irrespective  of  the  wavelet  used. 

^If  r,i  is  rearranged  according  to  projection  angles,  the  resulting  matrix  is  symmetric  block- Toeplitz. 
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4.1  Specialization  for  noisy  projection  data 


The  noisy  projection  data  is  represented  as 


y  =  (j)f_  +  n  (68) 

where  n  is  the  noise  vector  which  we  model  as  zero-mean,  white  and  Gaussian  with  variance  A'. 
Thus  n  ~  iV(0,  X'ImnxMn)-  This  results  in  the  input-output  equation  (compare  with  Equation  43) 


where 

=  PWn.  (70) 

Since  P  is  an  orthogonal  matrix  and  W~^  =  NW^,  u  N{0,  XImNxMn)  where 


(71) 


As  before,  we  can  partition  Equation  69  as  (compare  with  Equation  50) 


'  r,i 

rr^i 

4- 

iLs- 

r,2 

T,d 

L  i.w  J 

Esid) 

(72) 


Now,  by  invoking  the  earlier  assumptions  that  the  elements  of  the  matrix  r,2  are  mostly  negligible 
and  that  all  M  (DC)  elements  of  £(ci)  are  equal  if  the  object  is  entirely  covered  by  aU  M  projections, 
we  obtain  the  following  two  equations 


t.- 

fiM 

N 


-h 

cM^  +  i/Q 


(73) 

(74) 


where  ~  N{0,  MX).  From  Equation  74,  the  maximum  likelihood  estimate  of  c  is  n/{MN).  This 
is  the  same  as  the  value  of  c  we  had  used  earlier  in  the  case  of  no  noise  (Equation  88).  The  MAP 
estimate  of  from  Equation  73,  is  given  by 


|^_  =  argmax[P(^^_l3^_)]  (75) 

=  -f  (76) 

where  P(,-  is  the  prior  covariance  of  ~  iV(0,  P^,_).  We  as¬ 

sume  to  be  tmcorrelated  (or,  equivalently,  to  be  a  diagonal  matrix)  with  a  geometrically 
decreasing  variance  from  coarse  to  fine  scales.  Thus  the  diagonal  elements  of  P^,_  are  given  by 

mlJksn)  =  for  s  =  h,  h  +  1, . . . ,  /  (77) 

where  p  is  a  regularization  parameter.  This  implies  that  the  projection  data  mostly  influences  the 
reconstruction  of  coarse  scale  features  and  the  prior  model  influences  the  reconstruction  of  fine 
scale  features.  This  is  desired  since  the  high  frequencies  (i.e.  the  fine  scales)  are  mostly  corrupted 
by  noise.  Equation  76  can  now  be  simplified  to  (compare  with  Equation  60) 

Is-  =  ^ts-  +  ^  (78) 
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Figure  10:  Left:  A  32  X  32  Phantom.  Middle:  CBP  reconstruction  using  a  ramp  filter  and  M  = 
N  =  32.  Right:  Complete  multiscale  reconstruction. 


where 

A  =  (Tl,  + 

(79) 

and 

lo¬ 

ll 

ip 

(80) 

Thus  the  only  change  in  the  reconstruction  algorithm  from  the  noiseless  case  is  that  different  expres¬ 
sions  for  A  and  6,  given  by  Equations  79  and  80  respectively,  are  to  be  used  now  for  reconstructions 
from  noisy  data,  with  A  and  p  serving  as  regularization  parameters.  An  increased  regrdarization 
results  by  increasing  A  for  a  fixed  p,  or  by  increasing  p  for  a  fixed  A. 


5  Results 

5.1  Demonstration  of  the  reconstruction  method 


Figure  10  shows  a  phantom  along  with  the  CBP  and  the  complete  multiscale  reconstruction.  The 
size  of  the  phantom  is  32  X  32  and  the  imaging  geometry  is  defined  hy  M  =  N  =  32.  We  have 
used  a  ramp  filter  (i.e.  we  have  not  introduced  any  regidarization)  in  the  CBP  reconstruction. 
Thus,  for  this  complete  data  case,  the  CBP  and  the  multiscale  reconstructions  should  be  similar. 
This  is  precisely  what  is  seen  in  Fig  10  and  confirmed  in  Fig  11  which  shows  a  section  through 
the  reconstructions®.  It  is  observed  that  the  maximum  absolute  difference  in  the  reconstructed 
pixel  intensities  is  of  the  order  of  10“®  when  different  wavelets  are  used  to  reconstruct  this  same 
phantom. 

In  Figure  12  we  show  the  reconstruction  of  a  32  X  32  phantom  from  projection  data  collected  using 
M  =  N  =  S2.  The  D3  wavelet  is  used  for  mtdtiscale  reconstruction.  Figure  12  is  the  combined 
scale  reconstruction  where  the  reconstruction  at  a  particidar  scale  includes  the  contributions  from 

*The  discrepancy  between  the  CBP  and  the  multiscale  reconstructions  near  the  edges  is  due  to  the  fact  that  we 
are  reconstructing  on  a  circle  in  case  of  the  latter  rather  than  on  the  entire  square  field-of-view  as  in  the  former. 
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Figure  11:  A  horizontal  section  through  the  pheintom  (solid  line)  and  the  CBP  (circles)  and  mul¬ 
tiscale  (broken  line)  reconstructions. 


all  scales  coarser  to  and  including  that  scale.  In  Figure  13  (the  separate  scale  reconstruction)  we 
show  the  individual  contributions  at  different  scales. 


5.2  Reconstruction  of  fine  scale  features  by  approximating  r,i  as  an  identity 
matrix 

In  Figure  14,  we  show  the  CBP  reconstruction,  the  complete  multiscale  reconstruction  and  the 
finest  scale  contribution  from  the  reconstruction  of  a  32  x  32  phantom  using  the  Haar  basis.  We  use 
only  the  diagonal  elements  of  matrix  r,i  for  reconstruction  cind  assume  that  they  are  aU  equal  to 
1.  It  can  be  shown  with  relative  ease  that  the  complete  multiscale  reconstruction  in  this  case  is  the 
same  as  simple  back- projection  reconstruction  having  a  1/r  blurring.  We  can  see,  from  the  finest 
scale  contribution,  that  if  the  goal  is  edge  reconstruction,  it  is  enough  to  approximate  Fji  by  an 
identity  matrix.  This  reduces  the  computational  complexity  enormously.  The  edge  reconstruction, 
in  the  multiscale  framework,  only  involves  the  1-D  Haar  transform  of  the  strip  integral  data  and 
subsequent  back-projection  of  fine  scale  coefficients. 
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Figure  12:  Combined  sccde  reconstructions,  using  the  £>3  wavelet,  of  a  32  X  32  phantom  with 
M  =  iV  =  32.  Top  row,  left:  Phantom.  Top  row,  middle:  CBP  reconstruction  using  a  ramp  filter. 
Top  row,  right:  Multiscale  reconstruction.  Reconstruction  including  contributions  from  aU  scales 
at  and  coarser  to:  Middle  row,  left:  Scale  1.  Middle  row,  middle:  Scale  2.  Middle  row,  right:  Scale 
3.  Bottom  row,  left:  Scale  4.  Bottom  row,  middle:  Scale  5.  The  reconstructions  in  top  row,  right 
and  bottom  row,  middle  are  identical. 


5.3  Reconstruction  from  incomplete  data 

In  Fig  15  we  show  a  32  X  32  phantom  reconstructed  using  the  Haar  wavelet  and  M  =  5  and  N  =  32, 
i.e.  using  5  uniformly  spaced  angular  projections,  each  containing  32  parallel  strips.  This  is  the 
incomplete  data  case  and  hence  the  multiscale  reconstruction  is  expected  to  be  free  of  the  many 
artifacts  which  arise  in  the  CBP  reconstruction  due  to  the  space-invariance  assumption.  This  is 
exactly  what  is  seen  in  the  figure.  This  figure  adso  indicates  the  resolution-accuracy  tradeoff.  The 
reconstruction  at  Scale  3  (middle  row,  right,  in  the  figure)  does  not  enable  one  to  distinguish  between 
the  two  circles,  but  the  reconstruction  has  very  few  artifacts.  The  distinguishability  increases, 
accompanied  by  an  increase  in  the  artifacts,  at  Scede  5  (i.e.  full  multiscale  reconstruction).  Fig  16 
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Figitre  13:  Separate  scale  reconstructions,  using  the  wavelet,  of  a  32  X  32  phantom  with  M  = 
N  —  32.  Top  row,  left:  Phantom.  Top  row,  middle:  CBP  reconstruction  using  a  ramp  filter.  Top 
row,  right:  Multiscale  reconstruction.  Reconstruction  including  contributions  from  aU  scales  at  and 
coarser  to:  Middle  row,  left:  Scale  1.  Middle  row,  middle:  Sccde  2.  Middle  row,  right:  Scale  3. 
Bottom  row,  left:  Scale  4.  Bottom  row,  middle:  Scale  5. 


shows  another  excimple  of  the  incomplete  data  case.  The  parameters  here  are  the  same  as  in  Fig  15 
and  the  same  suppression  of  artifacts  is  seen. 

5.4  Reconstruction  from  noisy  data 

Figure  17  shows  a  phantom  along  with  the  unregularized  and  regularized  (MAP)  complete  mul¬ 
tiscale  reconstructions  after  adding  5  dB  noise  to  the  projection  data.  The  size  of  the  phantom 
is  32  X  32  and  the  imaging  geometry  is  defined  by  M  =  iV  =  32.  We  have  used  the  jDs  wavelet 
for  multiscale  reconstructions.  The  two  regularized  reconstructions  correspond  to  regularization 
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Figure  14:  From  left  (a)  32  x  32  Phantom,  (b)  CBP  reconstruction  with  M  =  iV  =  32  (c)  Complete 
multiscale  reconstruction,  and  (d)  The  finest  scale  contribution. 


parameters  A  =  1,  p  =  0.5  and  A  =  p  =  1  respectively.  Figure  18  shows  a  section  through  the 
reconstructions.  From  Figures  17  and  18  we  see  that,  as  is  expected,  an  increased  regtdarization 
results  when  the  value  of  p  is  increased  from  0.5  to  1,  keeping  A  constant  at  1. 


6  Conclusions 


We  have  used  the  natural  pixel  object  representation,  where  the  object  is  expanded  in  the  same 
basis  functions  along  which  the  strip  integral  data  is  collected,  as  the  starting  point  of  our  midtiscale 
reconstruction  method.  In  the  natural  pixel  representation,  the  natural  pixels  (i.e.  the  expansion 
coefficients)  and  the  strip  integral  data  are  related  by  the  system  matrix  G.  We  decompose  the 
natural  pixels  and  the  strip  integral  data  in  a  mtdtiscale  basis  so  that  they  are  related  by  a  new 
system  matrix  Fg  which  is  sparse  and,  in  addition,  has  a  symmetric  block-Toeplitz  structme  if 
the  M  angular  projections  are  uniform  between  0  and  180°.  Fast  inversion  adgorithms  exist  for 
these  matrices.  For  the  incomplete  data  case,  the  multiscale  reconstruction  is  relatively  free  of  the 
many  artifacts  that  plague  the  CBP  reconstruction.  We  have  shown,  through  examples,  that  for  an 
iU-posed  reconstruction  problem,  some  regularization  can  be  achieved  by  only  combining  coarser 
scale  reconstructions  instead  of  reconstructions  at  all  the  scales.  For  noisy  data  problems,  we  have 
specialized  our  multiscale  reconstruction  method  to  yield  MAP  multiscale  reconstruction  estimates 
corresponding  to  a  chosen  prior  on  the  multiscale  object  coefficients.  We  have  also  shown  that  a 
fast  method  of  detecting  fine  scale  features,  like  edges  and  boundaries,  in  an  object  is  by  taking  the 
Haar  transform  of  the  projection  data  and  subsequently  back-projecting  the  fine  scale  coefficients. 

Currently,  we  are  concentrating  on  regularizing  ill-posed  reconstruction  problems  by  incorporating 
different  stochastic  priors  in  our  multiscale  reconstruction  method.  We  use  the  multiscale  framework 
to  construct  regularizing  tree-based  models  in  scale  reflecting  stochastic  prior  information  about  the 
object.  Since  fast  algorithms  exist  for  performing  estimation  on  such  multiscale  trees  (i.e.  scale- 
recursive  estimation  algorithms  which  are  a  generadization  of  the  Rauch- Timg-Striebel  optimal 
smoothing  algorithm)  [14],  our  multiscale  reconstruction  method  enables  us  to  tackle  iU-posed 
reconstruction  problems  at  a  speed  much  faster  than  existing  methods.  In  addition,  our  technique 
provides  object  estimates  at  multiple  scales  along  with  corresponding  error  covariance  information. 
This  information  is  essentied  for  assessing  the  resolution-accuracy  tradeoff  and  determining  an 
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Figure  15:  Combined  scale  reconstructions,  using  the  Haar  wavelet,  of  a  32  x  32  phantom  with 
M  =  5  and  N  =  32.  Top  row,  left:  Phantom.  Top  row,  middle:  CBP  reconstruction  using  a  ramp 
filter.  Top  row,  right:  Multiscale  reconstruction.  Reconstruction  including  contributions  from  all 
scales  at  and  coarser  to:  Middle  row,  left:  Scale  1.  Middle  row,  middle:  Scale  2.  Middle  row,  right: 
Scale  3.  Bottom  row,  left:  Scale  4.  Bottom  row,  middle:  Scale  5.  The  reconstructions  in  top  row, 
right  and  bottom  row,  middle  are  identical. 


optimal  scale  for  reconstruction  of  each  portion  of  the  object. 
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Figture  16:  Combined  scale  reconstructions,  using  the  Haar  wavelet,  of  a  32  x  32  phantom  with 
M  =  5  and  N  =  32.  Top  row,  left:  Phantom.  Top  row,  middle:  CBP  reconstruction  using  a  ramp 
filter.  Top  row,  right:  Multiscale  reconstruction.  Reconstruction  including  contributions  from  all 
scales  at  and  coarser  to:  Middle  row,  left:  Scale  1.  Middle  row,  middle:  Scale  2.  Middle  row,  right: 
Scale  3.  Bottom  row,  left:  Scale  4.  Bottom  row,  middle:  Scale  5.  The  reconstructions  in  top  row, 
right  and  bottom  row,  middle  are  identical. 


Appendices 

Appendix  1;  The  matrix  r,d  is  circulant 

Pjd  has  the  form 


<  > 

<  ^ldl5JL2dl  > 

<  > 

<  ^2ciuJl2<il  > 

•••  <^2dlf^Mdl> 

<  ^MdU^ldl  > 

<  ^M<il)^2dl  > 

•••  <^Mdl^^Mdl> 
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Figure  17:  Prom  left  (a)  32  x  32  Phantom,  (b)  Unregularized  complete  multiscale  reconstruction 
with  M  =  N  =  Z2  and  5  dB  noise,  (c)  Regularized  complete  multiscale  reconstruction  with  A  =  1 
and  p  =  0.5,  and  (d)  Regularized  complete  multiscale  reconstruction  with  A  =  1  and  p  =  1. 


Figure  18:  A  horizontal  section  through  the  phantom  (solid  line)  and  the  various  reconstructions 
of  the  previous  figure. 
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where  <  e,g  >  refers  to  the  inner  product  or,  equivalently,  the  area  of  intersection  of  e  and  g. 
Clearly,  since  <  ^idukjdi  >=<  >,  is  synametric.  Also,  all  <  ^du^di  >■,  i  = 

1,2,  ...,M  are  equal  to  1  because  of  our  usage  of  normalized  basis  functions.  Further,  because 
of  our  assumption  of  uniformly  spaced  projection  angles,  <  ^di^^jdi  >  is  a  function  of  only 
|(i)mod(M)  —  (ji)mod(M)|.  This  results  in  a  symmetric  circulant 


Let  the  constant  row  sum  of  r,d  be  a.  Then 

M 

a  =  Yl<  ^idu^di  >  • 

t=i 

It  can  be  shown,  by  elementary  geometry,  that 


^  ^Idl  >  ^^dl  ^ 


sin  0  4-  cos  0  —  1 
sin  6  cos  6 

1 

<  ^IdU  ^M+2-.-)dl  > 


for  0  <  0  <  90°, 
for  e  =  0,  90°, 

for  i  =  1,2,...,  M/2  if  M  is  even,  and 
forz=  1,2,...,(M+ l)/2  if  Mis  odd. 


(82) 


(83) 


From  Equation  83,  it  can  be  shown  that  <  >  achieves  a  minimum  value  equal  to  2{\/^  - 

1)  =  0.83  when  i  corresponds  to  a  projection  at  45°  or  135°.  The  maximum  value  of  <  $idi,  $tdi  >  is 
1  when  i  corresponds  to  projections  at  0°  or  90°.  For  most  applications,  it  is  enough  to  approximate 
the  row  sum  of  F^d?  ci,  by  M.  For  M  =  16,  a  has  a  value  14.2,  and  for  M  =  32  a  value  of  28.3. 


Appendix  2:  Details  of  the  NP  based  multiscale  reconstruction  method 


By  adding  the  last  M  equations  in  the  matrix  partitioned  equation  (56),  we  get 


r.i 

V 

aM 

(84) 


Now,  by  applying  the  matrix  inversion  lemma  to  above  [13], 


£,-  =  Kit-  -  rji*!! 


-  if 
-  aM 


(85) 


Also,  from  Equation  84, 

^  =  (86) 

Now  the  elements  of  n  are  the  row  sums  of  Ffj.  Since  the  elements  of  F,2  are  the  areas  of  intersection 

of  basis  functions  k  =  1 . M}  with  k  =  s  =  h,  h  +  1, . . . n- 

1, . . . ,  n,},  they  are  mostly  negligible’^.  Thus  we  can  approximate 


^lM 

N 


K.  aMc 


(87) 


^Appendix  3  demonstrates  the  validity  of  this  assumption  for  the  Haar  case.  For  arbitrary  compactly  supported 
Daubechies  wavelets,  this  assumption  is  seen  to  be  valid  even  though  it  is  difficult  to  calculate  numerical  bounds  as 
is  done  in  the  Haar  case. 
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or 


H  ^  fi 
aN  ~  MN 


using  the  approximation  a  ^  M  (refer  to  Appendix  1). 


(88) 


Thus  the  multiscale  reconstruction  reduces  to  calculating  and  back-projecting 


i- 


(89) 


where  is  given  by  Equation  85.  Let  =  P#,  where  $  is  the  matrix  containing  the  multiscale 
basis  functions  arranged  according  to  projection  angles  eutid  P  is  the  same  permutation  matrix  as 
discussed  earlier.  The  multiscale  reconstruction  is  then  given  by 


+ 


MN 


sd±M- 


(90) 


Since  we  assume  that  the  object  is  fully  covered  by  all  M  projections,  we  are  only  interested  in 
reconstructing  the  intersection  of  the  field-of-views  of  the  M  projections.  Thus 


MN 


(91) 


Appendix  3;  The  elements  of  the  matrix  r,2  are  mostly  negligible  -  demonstration 
for  the  Haar  case 

The  elements  of  r,2  are  the  areas  of  intersection  of  basis  functions  k  =  1,. .  .,M}  with 

k  =  1,...,M-,  s  —  h,h+  n  =  1, . .  .,n,}.  The  maximum  absolute  value  in  r,2 

corresponds  to  <  >  (or  <  >),  where  projections  ki  and  k2  are 

separated  by  45°.  This  value  is  equal  to  5\/2/8  -  1  =  0.12.  But  the  majority  of  the  terms  in  r,2 
correspond  to  the  areas  of  intersection  with  fine  scale  basis  functions,  which  are  negligible.  As  an 
example,  for  fine  scales  s  =  h  and  h  -f  1,  the  absolute  value  of  <  >  is  bounded  by 


4  (sin  0  cos  0) 

where 

z  =  (0.5)^+^-* 

and  projections  i  and  j  are  separated  by  angle  9  with  0  <  0  <  90°. 

0  <  0  <  90°,  as  for  90°  <  9  <  180°  the  same  bounds  apply.  For  9 
intersection  <  >  are  identically  zero. 

Thus,  as  an  example,  for  9  =  11.25°,  >  has  a  maximum  absolute  value  of  0.02  for 

s  =  h  and  0.05  for  s  =  h  -|-  1.  For  9  =  45°,  these  values  are  0.008  and  0.02  respectively. 


(93) 

It  is  enough  to  consider 
=  0  or  90°,  the  areas  of 


1-t  1  + 


8z 


+ 


1622 


(sin^ -f  COS0  —  1)  (sin0 COS0  —  1) 


(92) 
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Appendix  4:  Justification  for  odd  M 


In  case  of  odd  M,  [rg2  r,d]  has  full  row  rank.  Now 


N-M  ~  [  ]  £/(d) 

(94) 

=  r,2i,_+r,4,(d). 

(95) 

Invoking  the  earlier  assumption  that  the  elements  of  r,2  are  mostly  negligible 

^Im  «  r,ai,(d) 

(96) 

or 

t(<!) 

(97) 

where  r  is  a  vector  containing  the  row  sums  of  Now,  as  proved  in  Appendix  5,  a  circulant 
matrix  and  it’s  inverse  both  have  constant  row  STuns.  Moreover,  the  row  sum  of  the  inverse  is  the 
reciprocal  of  the  row  smn  of  the  matrix  itself.  Since  is  circular  and  it’s  row  sum  a  «  M,  all  M 
elements  of  r  are  approximately  equed  to  1  /M.  Hence 

=  cljif.  (99) 

Thus  all  M  elements  of  i^{d)  are  equal  for  odd  M. 


Appendix  5 


Claim;  The  inverse  of  a  circulant  matrix  C,  if  it  exists,  has  constant  row  sums  which  are  reciprocal 
of  the  constant  row  sum  of  C. 


Proof;  Consider  a  M  X  M  circulant  matrix  C.  We  can  always  diagonalize  C  as  follows  [15] 

C  =  F*DF  (100) 

with 


F^F^  = 


y/M 


11...  1 

1  exp  [-j(27r/M)]  ...  exp  [-j(2x/M)(M  -  1)] 


(101) 


[l  exp  [— ji(27r/M)(M  -  1)]  ...  exp  [-y(27r/M)(M  -  1)^]  J 
where  *  denotes  complex  conjugation,  F*F  =  FF*  =  I  and  D  —  diag(<ii,  dj, . . . ,  iIm)-  Now 


CIm  =  F*DFli^f  =  F*  X  diag(di, d2,  ■  dja)  x 


y/M 

0 


(102) 
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since 

Af-l 

exp  [±j(2'K IM)ki]  =  0  for  ^  an  integer  and  1  <  ^  <  M  -  1. 

t=0 


Hence 


CIm  =  VMdiF* 


1 

0 


dikM- 


Now 

Hence 


d\ 


1, 


Thus,  from  Equations  104  and  106,  the  proof  is  complete. 


(103) 


(104) 


(105) 


(106) 


Appendix  6:  Summary  of  notations 

Tables  1  and  2  smnmarize  the  notations  developed  for  the  conventional  and  the  multiscale  recon¬ 
struction  techniques,  respectively. 
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Table  1:  Notations  developed  for  conventional  reconstruction  techniques. 


Notation 

Explanation 

M 

Number  of  angular  projections. 

N 

Niimber  of  strip  integrals  in  each  angular  projection. 

f 

f 

The  discretized  object,  defined  on  a  p  x  p  square  grid.  We  set  p  =  N. 

The  reconstructed  object. 

y{k) 

Projection  data  set  at  angle  ife,  ^  =  1, . . . ,  M. 

y{k)  =  [yki  yk2  •  •  •  ykNf- 

y 

Full  projection  data  set. 

y  =  y^(2)  •  • .  p^(M)]^. 

x{k) 

Object  coefficient  set  at  angle  k,  k  =  1,...,  M. 

X{k)  =  [Xkl  Xk2  ...  XkNf. 

X 

Full  object  coefficient  set. 

X  =  [®^(1)  ®^(2)  . . .  x^{M)]^. 

<!> 

The  projection  operator,  y  =  (l>f. 

The  back-projection  operator,  /  =  4>'^x. 

T 

-*-a 

The  matrix  representing  the  CBP  ramp  filtering  operation  for  a 
projection  data  set  at  a  particular  angle. 
x{k)  =  Tay{k),k  =  l,...,M. 

ImxM 

The  M  X  M  identity  matrix. 

0 

The  Kronecker  product. 

T 

T  =  ImxM  ®Ta,x  =  Ty. 

G 

The  NP  input-output  matrix,  y  =  Gx,  G  = 
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Table  2:  Notations  developed  for  mtdtiscale  reconstruction  technique. 


Notation 

Explanation 

ksn 

Multiscale  quantities  are  indexed  by  ksn. 

k  is  the  angle  of  projection,  k  —  1,...,  M. 

s  is  the  scale,  s  =  1, 2, . . .,  1  -f  IniV/ ln2 ;  iV  =  2*,i  G  Z'^. 

n  is  the  shift  within  the  scale,  n  =  1, 2, . . . ,  n* ; 

n,  =  iV/2'  for  s  =  1, 2, . . . ,  IniV/  In  2  and  n,  =  1  for  s  =  1  -|-  In N/  In 2. 

h 

The  finest  scjile,  s  =  1  =  h. 

1 

The  coarsest  scale,  s  =  hi  N/  In  2  =  /. 

d 

The  DC  component,  5  =  1  -f  lniV/ln2  =  d. 

Wa 

The  matrix  realization  of  the  1-D  wavelet  transform  operation. 

w 

W  =  ImxM  <S)  Wa. 

The  1-D  wavelet  transform  of  y{k),  }l^{k)  =  Way{k). 

ts 

%  =  Wy  = 

The  1-D  wavelet  transform  of  x{k),  ^(k)  =  Wax(k). 

is 

$ 

The  multiscede  projection  operator,  ^  =  #/,$  = 

The  multiscede  back-projection  operator,  £  = 

Te 

The  multiscale  input-output  matrix  arranged  according  to  projection 
angles.  ^  Fs  = 

P 

The  orthogonal  permutation  matrix  that  changes  projection  angle 
ordering  to  scale  ordering. 

3^(i) 

Contains  all  projection  data  at  scale  j. 
t^(j)  =  l^^(hj)^(2,j) 

Contains  projection  data  arranged  according  to  scales. 

it  =  +  1)  • '  •  3^ (0  ^ 

i.U) 

Contains  all  object  coefficients  at  scale  j. 

L 

Contains  object  coefficients  arranged  according  to  scales. 

is  =  ^is  =  f  (^ + 1)  •  •  •  f  (0 

r. 

The  multiscale  input -output  matrix  arranged  according  to  scales. 

±.  =  r,  =  PFeP^. 
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